home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / MLSS.FOR < prev    next >
Encoding:
Text File  |  1985-11-29  |  6.0 KB  |  200 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE MLSS
  5. C
  6. C        PURPOSE
  7. C           SUBROUTINE MLSS IS THE SECOND STEP IN THE PROCEDURE FOR
  8. C           CALCULATING THE LEAST SQUARES SOLUTION OF MINIMAL LENGTH
  9. C           OF A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH SYMMETRIC
  10. C           POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX.
  11. C
  12. C        USAGE
  13. C           CALL MLSS(A,N,IRANK,TRAC,INC,RHS,IER)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           A     - COEFFICIENT MATRIX IN FACTORED FORM AS GENERATED
  17. C                   BY SUBROUTINE MFSS FROM INITIALLY GIVEN SYMMETRIC
  18. C                   COEFFICIENT MATRIX A STORED IN N*(N+1)/2 LOCATIONS
  19. C                   A REMAINS UNCHANGED
  20. C           N     - DIMENSION OF COEFFICIENT MATRIX
  21. C           IRANK - RANK OF COEFFICIENT MATRIX, CALCULATED BY MEANS OF
  22. C                   SUBROUTINE MFSS
  23. C           TRAC  - VECTOR OF DIMENSION N CONTAINING THE
  24. C                   SUBSCRIPTS OF PIVOT ROWS AND COLUMNS, I.E. THE
  25. C                   PRODUCT REPRESENTATION IN TRANSPOSITIONS OF THE
  26. C                   PERMUTATION WHICH WAS APPLIED TO ROWS AND COLUMNS
  27. C                   OF A IN THE FACTORIZATION PROCESS
  28. C                   TRAC IS A RESULTANT ARRAY OF SUBROUTINE MFSS
  29. C           INC   - INPUT VARIABLE WHICH SHOULD CONTAIN THE VALUE ZERO
  30. C                   IF THE SYSTEM OF SIMULTANEOUS EQUATIONS IS KNOWN
  31. C                   TO BE COMPATIBLE AND A NONZERO VALUE OTHERWISE
  32. C           RHS   - VECTOR OF DIMENSION N CONTAINING THE RIGHT HAND SIDE
  33. C                   ON RETURN RHS CONTAINS THE MINIMAL LENGTH SOLUTION
  34. C           IER   - RESULTANT ERROR PARAMETER
  35. C                   IER = 0 MEANS NO ERRORS
  36. C                   IER =-1 MEANS N AND/OR IRANK IS NOT POSITIVE AND/OR
  37. C                           IRANK IS GREATER THAN N
  38. C                   IER = 1 MEANS THE FACTORIZATION CONTAINED IN A HAS
  39. C                           ZERO DIVISORS AND/OR TRAC CONTAINS
  40. C                           VALUES OUTSIDE THE FEASIBLE RANGE 1 UP TO N
  41. C
  42. C        REMARKS
  43. C           THE MINIMAL LENGTH SOLUTION IS PRODUCED IN THE STORAGE
  44. C           LOCATIONS OCCUPIED BY THE RIGHT HAND SIDE.
  45. C           SUBROUTINE MLSS DOES TAKE CARE OF THE PERMUTATION
  46. C           WHICH WAS APPLIED TO ROWS AND COLUMNS OF A.
  47. C           OPERATION IS BYPASSED IN CASE OF A NON POSITIVE VALUE
  48. C           OF IRANK
  49. C
  50. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  51. C           NONE
  52. C
  53. C        METHOD
  54. C           LET T, U, TU BE THE COMPONENTS OF THE FACTORIZATION OF A,
  55. C           AND LET THE RIGHT HAND SIDE BE PARTITIONED INTO A FIRST
  56. C           PART X1 OF DIMENSION IRANK AND A SECOND PART X2 OF DIMENSION
  57. C           N-IRANK. THEN THE FOLLOWING OPERATIONS ARE APPLIED IN
  58. C           SEQUENCE
  59. C           (1) INTERCHANGE RIGHT HAND SIDE
  60. C           (2) X1 = X1 + U * X2
  61. C           (3) X2 =-TRANSPOSE(U) * X1
  62. C           (4) X2 = INVERSE(TU) * INVERSE(TRANSPOSE(TU)) * X2
  63. C           (5) X1 = X1 + U * X2
  64. C           (6) X1 = INVERSE(T) * INVERSE(TRANSPOSE(T)) * X1
  65. C           (7) X2 =-TRANSPOSE(U) * X1
  66. C           (8) X2 = INVERSE(TU) * INVERSE(TRANSPOSE(TU)) * X2
  67. C           (9) X1 = X1 + U * X2
  68. C           (10)X2 = TRANSPOSE(U) * X1
  69. C           (11) REINTERCHANGE CALCULATED SOLUTION
  70. C           IF THE SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS IS SPECIFIED
  71. C           TO BE COMPATIBLE THEN STEPS (2), (3), (4) AND (5) ARE
  72. C           CANCELLED.
  73. C           IF THE COEFFICIENT MATRIX HAS RANK N, THEN THE ONLY STEPS
  74. C           PERFORMED ARE (1), (6) AND (11).
  75. C
  76. C     ..................................................................
  77. C
  78.       SUBROUTINE MLSS(A,N,IRANK,TRAC,INC,RHS,IER)
  79. C
  80. C
  81. C        DIMENSIONED DUMMY VARIABLES
  82.       DIMENSION A(1),TRAC(1),RHS(1)
  83.       DOUBLE PRECISION SUM
  84. C
  85. C        TEST OF SPECIFIED DIMENSIONS
  86.       IDEF=N-IRANK
  87.       IF(N)33,33,1
  88.     1 IF(IRANK)33,33,2
  89.     2 IF(IDEF)33,3,3
  90. C
  91. C        CALCULATE AUXILIARY VALUES
  92.     3 ITE=IRANK*(IRANK+1)/2
  93.       IX2=IRANK+1
  94.       NP1=N+1
  95.       IER=0
  96. C
  97. C        INTERCHANGE RIGHT HAND SIDE
  98.       JJ=1
  99.       II=1
  100.     4 DO 6 I=1,N
  101.       J=TRAC(II)
  102.       IF(J)31,31,5
  103.     5 HOLD=RHS(II)
  104.       RHS(II)=RHS(J)
  105.       RHS(J)=HOLD
  106.     6 II=II+JJ
  107.       IF(JJ)32,7,7
  108. C
  109. C        PERFORM STEP 2 IF NECESSARY
  110.     7 ISW=1
  111.       IF(INC*IDEF)8,28,8
  112. C
  113. C        CALCULATE X1 = X1 + U * X2
  114.     8 ISTA=ITE
  115.       DO 10 I=1,IRANK
  116.       ISTA=ISTA+1
  117.       JJ=ISTA
  118.       SUM=0.D0
  119.       DO 9 J=IX2,N
  120.       SUM=SUM+A(JJ)*RHS(J)
  121.     9 JJ=JJ+J
  122.    10 RHS(I)=RHS(I)+SUM
  123.       GOTO(11,28,11),ISW
  124. C
  125. C        CALCULATE X2 = TRANSPOSE(U) * X1
  126.    11 ISTA=ITE
  127.       DO 15 I=IX2,N
  128.       JJ=ISTA
  129.       SUM=0.D0
  130.       DO 12 J=1,IRANK
  131.       JJ=JJ+1
  132.    12 SUM=SUM+A(JJ)*RHS(J)
  133.       GOTO(13,13,14),ISW
  134.    13 SUM=-SUM
  135.    14 RHS(I)=SUM
  136.    15 ISTA=ISTA+I
  137.       GOTO(16,29,30),ISW
  138. C
  139. C        INITIALIZE STEP (4) OR STEP (8)
  140.    16 ISTA=IX2
  141.       IEND=N
  142.       JJ=ITE+ISTA
  143. C
  144. C        DIVISION OF X1 BY TRANSPOSE OF TRIANGULAR MATRIX
  145.    17 SUM=0.D0
  146.       DO 20 I=ISTA,IEND
  147.       IF(A(JJ))18,31,18
  148.    18 RHS(I)=(RHS(I)-SUM)/A(JJ)
  149.       IF(I-IEND)19,21,21
  150.    19 JJ=JJ+ISTA
  151.       SUM=0.D0
  152.       DO 20 J=ISTA,I
  153.       SUM=SUM+A(JJ)*RHS(J)
  154.    20 JJ=JJ+1
  155. C
  156. C        DIVISION OF X1 BY TRIANGULAR MATRIX
  157.    21 SUM=0.D0
  158.       II=IEND
  159.       DO 24 I=ISTA,IEND
  160.       RHS(II)=(RHS(II)-SUM)/A(JJ)
  161.       IF(II-ISTA)25,25,22
  162.    22 KK=JJ-1
  163.       SUM=0.D0
  164.       DO 23 J=II,IEND
  165.       SUM=SUM+A(KK)*RHS(J)
  166.    23 KK=KK+J
  167.       JJ=JJ-II
  168.    24 II=II-1
  169.    25 IF(IDEF)26,30,26
  170.    26 GOTO(27,11,8),ISW
  171. C
  172. C        PERFORM STEP (5)
  173.    27 ISW=2
  174.       GOTO 8
  175. C
  176. C        PERFORM STEP (6)
  177.    28 ISTA=1
  178.       IEND=IRANK
  179.       JJ=1
  180.       ISW=2
  181.       GOTO 17
  182. C
  183. C        PERFORM STEP (8)
  184.    29 ISW=3
  185.       GOTO 16
  186. C
  187. C        REINTERCHANGE CALCULATED SOLUTION
  188.    30 II=N
  189.       JJ=-1
  190.       GOTO 4
  191. C
  192. C        ERROR RETURN IN CASE OF ZERO DIVISOR
  193.    31 IER=1
  194.    32 RETURN
  195. C
  196. C        ERROR RETURN IN CASE OF ILLEGAL DIMENSION
  197.    33 IER=-1
  198.       RETURN
  199.       END
  200.